Assessing the potential of composite confining systems for secure and long-term CO2 retention in geosequestration

A potential geologic target for CO2 storage should ensure secure containment of injected CO2. Traditionally, this objective has been achieved by targeting reservoirs with overlying seals-regionally extensive, low permeability units that have been proven capable of retaining buoyant fluid accumulations over geologic time. However, considering that the amount of CO2 is limited by a decadal injection period, vertical migration of CO2 can be effectively halted by a composite system of discontinuous shale/silt/mudstone barriers in bedded sedimentary rocks. Here, we studied the impact of depositional architectures in a composite confining system on CO2 migration and confinement at reservoir scale. We stochastically generated lithologically heterogeneous reservoir models containing discontinuous barriers consistent with statistical distributions of net-sand-to-gross-shale ratio (NTG) and horizontal correlation lengths derived from well log data and observations of producing hydrocarbon fields in Southern Louisiana. We then performed an extensive suite of reservoir simulations of CO2 injection and post-injection to evaluate the sensitivity of CO2 plume migration and pressure response of the composite system to a series of geologic and fluid parameters including the lateral continuity of barriers, NTG, permeability anisotropy within the sand body, and capillary pressure contrast between the sand and shale facies. The results indicate that lateral continuity of barriers and NTG are the dominant controls on CO2 plume geometry and pressure build-up in the reservoir, while the impact of NTG is particularly pronounced. The significance of intraformational barriers becomes apparent as they facilitate the local capillary trapping of CO2. Those barriers improve the pore space occupancy by promoting a more dispersed shape of the plume and ultimately retard the buoyancy-driven upward migration of the plume post injection.

beneath the flow barriers is called capillary heterogeneity trapping induced by the capillary pressure contrast between the sandbody and fine-grained barriers 20,21 .Local trapping of CO 2 under the intraformational barriers is equivalent to the capillary pinning mechanism 22,23 that retards the buoyancy-driven migration of CO 2 plume to shallow depths of the reservoir and with a good design, the composite confining system will prevent CO 2 from reaching the top of the confining system.In addition, the tortuous flow path of CO 2 plume created by the existence of discontinuous barriers causes a large volume of rock and resident brine to be in contact with CO 2 , and hence, facilitates other CO 2 trapping mechanisms such as snap-off, dissolution, and mineral trapping.Capillary heterogeneity trapping of CO 2 controlled by intraformational barriers has been demonstrated as an effective mechanism for CO 2 containment in previous studies 15,17,18,24 .
The heterogeneity in the spatial distribution and lateral continuity of discontinuous shale barriers as well as sand/shale ratios are among the major uncertainty associated with the composite confining systems 5,25 .These stratigraphic heterogeneities cannot be resolved deterministically using seismic data.Additionally, quantifying the lateral heterogeneity of shale barriers using well-to-well correlation data is challenging as the data are sparse due to the large spacing between the wells.Thus, correlating the stratigraphic distribution of barriers between the wells using deterministic modeling is difficult and interpretations are not typically distinctive 25 .Efforts have been undertaken to utilize stochastic techniques in lieu of conventional mapping to model complex heterogeneity distribution within a heterogeneous reservoir 25,26 .The main intention of this study is to consider a wide range of lateral continuity and density of shale barriers in geologic models and investigate the impact of these architectures and their characteristics on the migration and confinement of CO 2 within these formations.To this end, we employed a stochastic modeling approach to generate heterogeneous facies models containing horizontal discontinuous shale barriers populated in sand bodies, controlled by the correlation length of barriers and sand/shale ratios as factors constrained by Spontaneous Potential (SP) well-log data and observations of producing hydrocarbon fields 5 .
In this study, we performed a set of reservoir simulations involving CO 2 injection and post-injection in a composite confining system to evaluate the sensitivity of CO 2 plume migration and pressure response of the system to a range of geologic and fluid parameters.Net-sand-to-gross-shale ratio (NTG) and horizontal correlation lengths of shale barriers were considered as factors controlling the degree of heterogeneity in the geologic models.The composite confining models generated in this study serve as representative examples of Southern Louisiana Miocene deltaic deposits with large numbers of interbedded barriers and no regionally extensive caprock seals.To simplify the model, we did not account for a distinct injection zone commonly considered in actual storage projects.Instead, we injected CO 2 into the lower section of the confining zone.Using the reservoir simulation results, the performance of a composite system for CO 2 confinement has been evaluated and the corresponding sensitivity analysis results have been presented in the following sections.

Geostatistical model
We used a synthetic geologic model to perform a sensitivity analysis for characterizing CO 2 migration and confinement in a composite confining system.To mimic a composite system of flow barriers, we have geostatistically generated two-dimensional (2D) heterogeneous facies models containing horizontal discontinuous shale barriers (i.e., capillary barriers) populated in sand bodies.The models involve 150 × 150 grid blocks with horizontal and vertical resolutions of 219 ft and 2.2 ft, respectively.The shale barriers were considered to be laterally correlated with correlation lengths ranging from 1973 ft (600 m, equivalent to 6% of the model's horizontal extent) to 26,212 ft (8 km, equivalent to 80% of the model's horizontal extent).The stochastic facies models were constrained by two metrics including statistical distribution of net-to-gross ratio (NTG) and horizontal correlation length of shale barriers.The distribution of NTG ratio was derived from evaluation of Spontaneous Potential (SP) logs from 323 wells in Southern Louisiana that record over 40,000 individual mudstones within the Miocene deltaic and shore zone facies.Log values were normalized to a range of −80 to +20 MV and a cutoff value of − 30 MV was employed to distinguish between sandy and muddy facies.The footprint of the hydrocarbon-water contacts was used to estimate the range of barrier lengths 5 .
The spatial distribution of shale facies within the 2D models is generated using a spectral method 27 .First, we generated a white random noise matrix of R( − → r ) with the size of 150 × 150.Then, we calculated its Fourier transform F ( − → r ) , multiplied by a Gaussian function: By taking the inverse Fourier transform of G ( − → s ) , we obtained a correlated Gaussian distribution field, R ( − → r ) , yielding an auto-correlation function: Where * denotes the convolution product and = π/s 0 represents the correlation length considered for the shale facies ( cl x ).For a given shale/sand ratio (NTG), the grid blocks corresponding to shale facies are defined as R ( − → r The horizontal permeability field for each facies was generated using a log-normal distribution using the mean values of 500 mD and 0.01 mD for sand and shale, respectively.The standard deviation of the natural log transform of permeability was assumed to be 0.45 for both facies.The permeability field of sand was considered anisotropic by taking the vertical-to-horizontal permeability ratio ( K v /K h , where K v and K h are the vertical and (1) horizontal permeability, respectively) as one of the uncertain parameters in this study.The porosity ( φ ) distribu- tion is assumed to be dependent on permeability as follows 28 : The spatial distribution of porosity and permeability in two instances of stochastic geologic models with different NTG and cl x is displayed in Fig. 1.

Reservoir model
We used an Equation-of-State (EoS) reservoir simulator for compositional modeling, namely, GEM TM [Computer Modelling Group Ltd. (CMG), Calgary, Alberta, Canada] 29 to simulate CO 2 storage in saline aquifers.Twodimensional reservoir models with dimensions of 32,850 ft and 330 ft in the horizontal and vertical directions are generated using a stochastic method.The reservoir's lateral dimension is taken to be sufficiently large for better monitoring of the CO 2 plume lateral extension in a long-term simulation.The upper and lower boundaries of the reservoir were considered closed, while open boundaries are imposed on the lateral sides of the model by extending the pore volume of the outer boundary grid blocks.The model has a dip of 1 • and depth of 8000 ft.The initial pressure and temperature of the reservoir were assigned using a hydrostatic gradient of 0.465 psi/ft and a geothermal gradient of 0.0165 F/ft, with surface pressure and temperature set at 14.7 psi and 60 F, respectively 30 .
The Brooks-Corey model was used to generate the relative permeability curves of water and CO 2 ( k rw , k rg ), given by 31,32 (3) K h = 7 × 10 7 (φ 9.61 ) , (4) Where S w denotes the water saturation.The irreducible water saturation ( S wi ) and the critical CO 2 saturation ( S gr ) are taken to be 0.2 and 0.05, respectively.To account for CO 2 residual trapping, we incorporated hysteresis characteristics into the relative permeability curve.To model the gas relative permeability hysteresis, we applied the Land trapping model, which relates the trapped gas saturation S gt to the initial gas saturation S gi during the imbibition as 33 Where C is the Land's trapping coefficient calculated as follows: Where S max g is the maximum gas saturation associated with the imbibition curve, and S max gr is the maximum residual gas saturation.In this study, the S max gr is assumed to be 0.2, and S max g = 1 − S wi .We did not quantify the uncertainty surrounding the impact of relative permeability on CO 2 plume geometry and reservoir pressure.Nevertheless, previous studies have demonstrated the effect of relative permeability parameters on pressure accumulation and CO 2 plume configuration in homogeneous reservoirs 34,35 .
The drainage capillary pressure was created using the Brooks-Corey model as follows: Where P ce is the capillary entry pressure.P ce is considered as one of the uncertain parameters, ranging from 0.03 psi to 60.2 psi for sand ( P ce,sand ) and 137 psi to 2146 psi for shale ( P ce,shale ) 36 .We did not account for hysteresis in the capillary pressure curves.
To consider the effect of local heterogeneity on capillary pressure, we employed the Leverett J-function 37 Where σ is interfacial tension and θ is contact angle.Using the J-function, we scaled the capillary pressure curve for each grid block using their porosity and permeability values, while σ and θ are considered spatially invariant.

Design of experiments
We conducted CO 2 injection simulations in geologic models by considering 100 different combinations of uncertain parameters, namely, NTG, cl x , K v /K h , P ce,sand , and P ce,shale .The range of variations for each uncertain param- eter is shown in Fig. 2. For all simulations, supercritical CO 2 was injected through a single well for a period of 10 years with an injection rate of 0.7 mmscf/day and a maximum well's bottom-hole pressure of 6400 psi.The maximum allowable injection pressure was selected to be 80% of the rock fracture pressure, which is estimated to be 8000 psi based on the lithostatic pressure gradient of 1 psi/ft.The injection rate and period were chosen to manage maximum pressure build-up and maintain stabilized injectivity.The injection well is located in the center of the domain with a perforation length of 44 ft equivalent to grid blocks of 130 to 149 in the vertical direction.
To simplify the model, we did not account for a distinct injection zone commonly considered in actual storage projects.Instead, we injected CO 2 into the lower section of the confining zone.The simulations continued for 200 years after injection stopped.One criterion for selection of 200-year post-injection simulations was to ensure the stabilization of CO 2 upward migration.We quantified the sensitivity of CO 2 plume migration and pressure response of the reservoir to fluid/geologic input features.The vertical movement of the CO 2 plume is quantified with the height of the center of mass of the plume (in the y direction) by calculating the first spatial moment, defined as 38 Where ρ is the CO 2 density and s denotes CO 2 saturation in each grid block.x and y are the coordinates of grid blocks and L and H are respectively the lateral extent and thickness of the geologic model.
The response metrics were selected as the average grid-based saturation of CO 2 plume ( S CO 2 ), the maximum dimensionless lateral ( L x /L ) and vertical ( L y /H ) extent of the CO 2 plume, where L x and L y are the maximum extent of the plume in the x and y directions, respectively, dimensionless center of mass of the CO 2 plume (6)

Results and discussion
CO 2 saturation and pressure plume dynamics We investigated the migration of CO 2 plume and pressure response of the confining system under various scenarios considering different combinations of uncertain variables.Through this parametric sensitivity analysis, we were able to determine the combined role of uncertain variables on CO 2 saturation and pressure plume behavior.Figure 3 displays CO 2 saturation and pressure profile at the end of injection and 200 years of post-injection for two geologic models with various NTG and cl x values, representing how heterogeneity can affect the shape of the plume and distribution of pressure buildup in the reservoir.Figure 4 represents the probability distribution of response variables including S CO 2 , L y /H , L x /L and CM/H obtained from injection and post-injection simulations using 100 combinations of uncertain parameters.The average CO 2 saturation of the plume was found to be 0.34 and 0.24 at the end of the injection and post-injection (200 years) stage, respectively.According to Fig. 4a, the distribution of S CO 2 corresponding to the post-injection has shifted towards lower values and has a long skew towards higher values, indicating the spread or redistribution of CO 2 plume and occupying a larger pore volume of the pores.As shown in Fig. 3, CO 2 plume tends to become channelized underneath the capillary barriers and locally trapped with a saturation higher than the CO 2 residual saturation.These large saturation regions of CO 2 structurally trapped underneath the barriers by capillary heterogeneity potentially contain mobile CO 2 17 .This type of trapping is as effective as structural trapping unless the capillary integrity of those barriers is compromised through a large aperture feature such as an open wellbore.
Figure 4b and c display the probability of the lateral and vertical migration of CO 2 plume at the end of injection and late post-injection.According to those graphs, we observed that CO 2 plume propagates both laterally and vertically during the post-injection stage.However, the variation in its lateral expansion is more pronounced than in the vertical direction.This means that while the buoyancy effects promote the vertical migration of CO 2 plume at post-injection, the existence of barriers exerts a significant constraint on the influence of buoyancy forces and hence the vertical growth of the plume.In addition, according to Fig. 4b, there were a few case scenarios, in which the CO 2 plume reached the top of the formation ( L y /H = 1).Increasing the thickness of the confining zone would likely increase the confinement capacity of those failed scenarios.In most scenarios, the plume mainly accumulated in the deeper parts of the confining zone, as the distributions have a long tail towards higher L y /H values.As shown in Fig. 4d, the location of the center of mass of the plume (CM) for all realizations did not exceed 50% of the reservoir thickness, meaning that CO 2 saturation tends to decrease towards the upper layers of the formation and the majority of the plume mass is concentrated in the deeper parts.Thus, the barriers facilitate the local trapping of CO 2 and retards its buoyancy-driven upward migration at the late post-injection.We also explored the correlation between the plume height and its lateral extension and showed how the center of mass of the plume varies with these two parameters (see Fig. 5).A negative correlation between the plume height ( L y /H ) and its lateral extent ( L x /L ) and positive correlation between L y /H and CM/H were observed.
In addition to the analysis of the CO 2 saturation plume (i.e., the plume configuration), we further investigated pressure propagation in the reservoir at the end of injection and 200 years post-injection.Figure 6 represents the probability distribution of maximum pressure buildup ( P max ), the pressure buildup at the top of the res- ervoir ( P top ), and the location of the pressure buildup occurrence ( H pmax /H , where H pmax is the distance of the location of the maximum pressure buildup occurrence from the bottom of the reservoir) obtained from the injection and post-injection simulations using 100 realizations of uncertain parameters.As can be seen in Fig. 6a, P max distributions show the same trend for both end of injection and post-injection, as they both skew toward larger values.The average pressure buildup at the end of injection and post-injection was 500 psi and 228 psi, respectively.Additionally, according to the distributions shown in Fig. 6a, pressure attenuates significantly at the late post-injection, manifesting approaching an equilibrium state for the reservoir's pressure after injection stops.As shown in Fig. 6b, the pressure buildup at the top of the reservoir is notably smaller than the maximum pressure buildup in the reservoir, meaning that the existence of intraformational barriers suppresses the vertical communication of the pressure in the reservoir.This can be better understood from Fig. 6c, showing that the maximum pressure buildup was observed at locations far from the top of the reservoir.Furthermore, we observed a sharp attenuation of P top after injection stops (see Fig. 6b). Figure 7 illustrates the variation of pressure along the reservoir depth at the end of injection and post-injection for one of the scenarios.It can be noted that the pressure plume is mainly accumulated at the greater reservoir depths, while the pressure of the upper layers remains almost unchanged during CO 2 storage.

Sensitivity analysis
We also examined the level of the sensitivity of CO 2 plume shape and pressure response of the confining zone to uncertain parameters including NTG, cl x , K v /K h , P ce,sand , and P ce,shale .A visual display of the sensitivity analysis is depicted in the tornado chart in Fig. 8.This chart represents the contribution of uncertain parameters to the response variables including S CO 2 , L y /H , L x /L , CM/H, P max , and P top .To build the tornado chart, one parameter at a time was varied between its lower and upper values taken from their uniform distribution, while other variables were kept at their median values.Then, the model response were laid out in the tornado chart.The degree of sensitivity of the uncertainty metrics to the input variables can be quantified by the bar lengths corresponding to each variable.As can be seen, the permeability anisotropy ( K v /K h ) of sand has no effects on the vertical migration of CO 2 plume ( L y /H and CM/H), which is yet highly controlled by the correlation length ( cl x ) and density (NTG) of capillary barriers.Although buoyancy-driven migration of CO 2 in homogeneous reservoir models has been shown to be highly influenced by the permeability anisotropy in previous studies 39 , our results delineate that the lateral extension and density of capillary barriers have more influential impact on retarding the upward migration of CO 2 .According to the tornado chart, it was observed that P ce,sand has minimal impacts on the configuration of CO 2 plume (i.e., L y /H , L x /L , and CM/H) and pressure buidlup (i.e., P max , and P top ).This finding suggests that the capillary entry pressure of sand bodies does not play a key role in controlling the migration of CO 2 plume and pressure propagation in a heterogeneous reservoir with capillary barriers.However, it is an influential parameter controlling CO 2 plume geometry in homogeneous reservoirs 40 .It was observed that lower values of P ce,shale promote the upward migration of CO 2 plume (i.e., L y /H ), leading to reduced lateral migration (i.e., L x /L ).In general, according to the sensitivity analysis, it is evident that the CO 2 plume geometry and pressure propagation are highly influenced by NTG and cl x in a heterogeneous medium.The box plots shown in Fig. 9 represent the dependence of CO 2 saturation ( S CO 2 ) and aspect ratio of the plume ( AR plume , which is L x L y ), on 1-NTG and cl x /L .According to Fig. 9a, we observed a positive correlation between the median value of S CO 2 with the length of barriers and their volume fraction in the confining zone.However, S CO 2 dependence on NTG appears to be more pronounced than that on cl x /L .According to Fig. 9c and d, AR plume median values show a positive correlation with NTG and cl x /L , while the correlation with NTG is more signifi- cant.The findings indicate that there are no substantial changes in the aspect ratio and spatial distribution of the plume when the lateral extent of the barriers exceeds 40% of the formation's length.
The box plots shown in Fig. 10 represent the dependence of P max on 1-NTG and cl x /L at the end of injection and 200 years post-injection.We observed an increasing trend in the maximum pressure buildup as the barrier length and their volume fraction increase, while the correlation with NTG is more pronounced.Considering the sensitivity of the formation's response to NTG, we recommend that the characterization of a composite confining system of heterogeneous barriers should be heavily focused on the quantification of this characteristic.
Figure 11 represents the variations in the plume aspect ratio ( AR plume ) as a function of maximum pressure buildup P max as well as the role of NTG and cl x /L on their correlation.According to these results, high aspect ratios correspond to high pressure buildups in the formation and they appear at large barrier lengths and large volume fractions of barriers.Although large barriers or a high density of barriers can promote the lateral extension of the plume, they can also lead to a higher pressure buildup, as pressure communication in the vertical direction is limited in the presence of barriers.However, according to Fig. 6, it should be noted that pressure accumulates mainly in the deeper parts of the confining zone.

Summary
This study investigates the containment of CO 2 in regards to its vertical migration within lithologically heterogeneous reservoir rocks composed of stacked fine-grained barriers interbedded within sand units known as composite confining systems 5 .The reservoir-scale two-phase flow simulations presented here are representative  of the pressure-driven and buoyancy-driven flow of CO 2 occurring during injection and post-injection in saline aquifers.Heterogeneous facies models were generated through geostatistically populating horizontal discontinuous shale barriers into sand bodies.Numerous geologic models were generated encompassing various combinations of uncertain parameters, including NTG, cl x , K v /K h , P ce,sand , and P ce,shale .We studied the impact of heterogeneous composite systems on the pressure-and buoyancy-driven migration of CO 2 , the effectiveness of confinement, and the formation's pressure response.
According to the simulation results, we observed that the existence of barrier layers exerts a significant influence on the vertical migration of the plume and the pressure propagation in the composite system.The dispersed flow path of CO 2 plume, controlled by the spatial heterogeneity of capillary barriers, is more pronounced after the cessation of injection when the plume migrates upward due to the effect of buoyancy forces.Based on the parametric sensitivity analysis, it has been demonstrated that the lateral continuity of discontinuous barriers (i.e., cl x ) and their volume fraction (i.e., NTG) play significant roles in determining the shape of the CO 2 plume and vertical pressure communication in the system.With increasing the fraction of barriers and their horizontal length, the vertical movement of CO 2 plume becomes restricted while its lateral movement is encouraged.Furthermore, the plume shape and pressure response and overall retention capacity of a composite system were found to be more susceptible to changes in NTG than cl x .Hence, NTG should be prioritized as the key factor for geologic characterization of composite systems, given its significant influence on the dynamic response of such systems to CO 2 injection.
With regard to the findings of this study, it is reasonable to suggest that formations containing capillary barriers are effective in containing the injected CO 2 plume within the formation, as heterogeneity serves to limit the reliance of the formation seal as the only mechanism for containment.Even though the composite system considered in this study acts as the injection zone, the efficacy of such a system for CO 2 containment at the post-injection stage, where plume buoyant migration is dominant, has been demonstrated.A study by Bump et al. 5 has provided further demonstration of the effectiveness of composite systems as a confining zone.The 3D reservoir-scale flow simulation in their study has shown that the migration of CO 2 plume could be predominantly horizontal with minimal intrusion into the confining zone with the extent of the capillary barriers occupying at least 6 % of the reservoir domain.

Figure 1 .
Figure 1.Two examples of geostatistical models with different net-to-gross (NTG) ratios and correlation lengths of impermeable barriers in the horizontal direction ( cl x ).The model extends 330 ft and 32,850 ft in the vertical and horizontal directions, respectively.The upper figures represent the facies distributions.Blue is sand and red represents impermeable barriers.The middle figures represent the porosity distributions.The lower figures show permeability distributions.
, y)ρ(x, y)s(x, y)x i y j dx dy (13) CM = M 01 M 00 Vol.:(0123456789) Scientific Reports | (2023) 13:21022 | https://doi.org/10.1038/s41598-023-47481-2www.nature.com/scientificreports/(CM/H), maximum pressure buildup in the reservoir ( P max ), the pressure buildup at the reservoir top ( P top ), and the dimensionless depth at which the maximum pressure buildup occurs ( H max /H ).The response variables were measured at the end of injection and 200-year post-injection.Our intention was to capture the viscous-and buoyancy-driven migration of CO 2 occurring during injection and post-injection, respectively.Thus, the simulations were conducted for 10 years of injection and continued for 200 years after injection stopped.

Figure 2 .
Figure 2. Paired scatter plot illustrating the distribution of five uncertain parameters including NTG, cl x , K v /K h , P ce,sand , and P ce,shale , used for the sensitivity analysis.The x-axis and y-axis of each plot represent two different variables of interest.100 different combinations of these uncertain parameters were sampled from their distributions to create various stochastic reservoir models.

Figure 3 .
Figure 3. S CO 2 saturation and pressure profile at the end of injection and 200 years of post-injection for two geologic models shown in Fig. 1.The figures in the left box (a-d) correspond to NTG = 0.82 and cl x = 3379 ft and the figures in the right box (e-h) correspond to NTG = 0.61 and cl x = 17,898 ft.L x and L y are the maximum extension of S CO 2 plume in the horizontal and vertical directions, respectively.

Figure 4 .
Figure 4. Probability distribution of output variables including the average CO 2 saturation ( S CO 2 ), dimensionless height of the plume ( L y /H ), dimensionless lateral extension of the plume ( L x /L ), and center of mass of the plume (CM/H) at the end of injection and 200 years post-injection.

Figure 5 .
Figure 5. Correlation between the plume height ( L y /H ) and its lateral extension ( L x /L ) and their relation to the plume center of mass (CM/H).

Figure 6 .
Figure 6.Probability distribution of the pressure-related output variables including the maximum pressure buildup ( P max ), the pressure buildup at the top the composite system ( P top ), and the location of the pressure buildup occurrence ( H pmax /H ) at the end of injection and 200 years post-injection.

Figure 7 .
Figure 7. Pressure profile along the depth of the composite system in one of the case scenarios at the end of injection and 200 years post-injection.

Figure 9 .
Figure 9. Box plots representing the dependence of CO 2 saturation ( S CO 2 ) and aspect ratio of the plume ( AR plume : which is L x L y ) on NTG and cl x /L.

Figure 10 .
Figure 10.Box plots representing the dependence of P max on NTG and cl x /L at the end of injection and 200 years post-injection.

Figure 11 .
Figure 11.Variations in the aspect ratio of the plume ( AR plume ) as a function of maximum pressure buildup ( P max ) as well as NTG and cl x /L.